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A spherical blast wave with relativistic velocity can be described by a similarity solution, 
that is used for theoretical models of gamma-ray bursts. We consider the linear stability 
of such a relativistic blast wave propagating into a medium with density gradient. The 

^ | perturbation can also be expressed by a self-similar form. We show that the shock front is 

^ , unstable in general, and we evaluate the growth rate. 

*0 ' 81. Introduction 

Gamma-ray bursts have puzzled astronomers since their accidental discovery in 
the late sixties. Recent observations still provide new topics of study, afterglow, 
optical flash, & brightest events, & and so forth. In these bursts, huge amounts 
of energy are released instantaneously. The subsequent evolution of the radiation- 
dominated matter can be described by an expanding blast wave. Gamma-rays are 
emitted when a matter with highly relativistic velocity, with Lorentz factor r > 10 2 



is converted into radiation. 



The observed fluctuations in gamma-rays with a short timescale are thought to 



be associated with the complex behavior of shock waves. Two locations have been 



proposed as the origin of the irregularities. One location is inside the relativistic 



fluid, as described by the internal shock model. When shock waves intersect each 
other, gamma-rays are emitted. The complexity of the corresponding time profile 
depends on the internal shock structure. The other location is in the shock front 
itself, as described by the external shock model. When the shock front impinges upon 
^ ■ ISM (inter stellar matter), gamma-rays are emitted. Both models have advantages 

and disadvantages. Waxmani & showed synchrotron emission from external shocks 
provides a successful model for the broken power law spectra and smooth temporal 
behavior of afterglow. Sari and Piran,l9 however, showed that external shocks cannot 
produce variable gamma-ray bursts unless they are produced by extremely narrow 
jets or if only a small fraction of the shell emits radiation. The expanding shell 
and ISM were assumed to be smooth. The observed nature with rapid variability 
depends on the fluctuations of the external matter, which is not yet known. The 
possibility of irregularity in the shock front can be investigated. An irregular shell 
with angular fluctuations in the external shock may produce temporal behavior of 
the gamma-ray emission. 

A theoretical model of the blast wave with ultra-relativistic speed was considered 
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by Blandford and McKee.i) They obtained a spherical self-similar solution, which 
can be regarded as the relativistic version of the Sedov- Taylor solution. Using a 
numerical code for the spherical symmetry, Kobayashi, Piran and SaritP simulated an 
expanding shock with relativistic speed and confirmed that the self-similar solution is 
a good approximation. The spherical symmetry may be too great a simplification in 
some cases. Indeed, in non-relativistic hydrodynamics, Ryu and Vishniac™ showed 
that under certain conditions the spherical shock becomes overstable with respect to 
perturbations with tangential velocity. From this point of view, it is important to 
examine the dynamical stability of the Blandford-McKee solution. 

In this paper, we consider the linear stability with respect to non-spherical de- 
formations of ultra-relativistic shock waves. There is a similarity solution describing 
shocks with an impulsive energy supply. In Section 2, we show that the perturbation 
functions can be obtained in a power low form. Matching the boundary conditions at 
the shock front and using a regularity condition at the origin, the temporal behavior, 
i.e. an eigenvalue of the functions, is determined. Section 3 is devoted to discussion. 

§2. Basic equations and approximate solutions 

2.1. Sphrical shock waves 

Relativistic fluid motion is described by conservation laws of number and energy- 
momentum. We explicitly give these equations in spherical coordinates for two- 
dimensional flow with the velocity (y r ,vg) and the Lorentz factor 7 as 

d t D + \d r {r 2 Dv r ) + —^— -d e (sm9Dv e ) = 0, (24a) 
r z r sin ft 

d t (W -p) + \d r (r 2 Wv r ) + —^—d e (sm9Wvg) = 0, (2-lb) 

dt (W v r ) + \d r {r 2 (W v 2 + p) } + — |— d g (sin 9 W v r vg) 
r A r smU 

v 2 --^ = 0, 2-lc 

r r 

8t(W Vg) + - d r (r W V r Vg) + — |— Og {sin 9 (W Vq + p) } 
r r smd 



+-Wv r vg-- cat = 0, (2-ld) 



where 



D = n 7 , (2-2) 
W = (e+p)j 2 = 4p 7 2 , (2-3) 

and n, e and p are the number density, energy density and pressure. We have assumed 
that after a strong shock, that the relativistic fluid is radiation dominated, i.e. e = 

Blandford and McKee& obtained a spherically symmetric similarity solution 
of Eqs. (2- la) - (2-ld). This solution can be written in terms of t and r. The 
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kinematical energy of the shock decreases with time. The Lorentz factor r of the 
shock front therefore indicates the timescale. In place of the radial coordinate, they 
introduced a new coordinate x as the similarity variable. The position of the shock 
front corresponds to x = 1, and the interior region corresponds to x > 1- These new 
coordinates are mathematically defined as 

r 2 = (t/t y m , (2-4) 

X = {l + 2(m + l)r 2 }(l-r/t), (2-5) 

where to is a constant, and the index m determines the energy supply rate of the blast 
wave. For the m = 3 case, the motion corresponds to an adiabatic point explosion, 
that is, an impulsive injection of energy at the center. For the case m < 3, the 
decrease of the velocity is less than that for the adiabatic case. This means that the 
solution corresponds to a blast wave with additional power supply. B This solution 
has additional inner shock and a contact discontinuity inside the flow. 

We assume a strong shock, that is, that the kinematical shock energy is extremely 
large. This situation corresponds to a large value of r. The solution matched with 
the strong shock condition can be expressed by 

Do = 2 ni r 2 h , (2-6) 
Po = \w x r 2 f , (2-7) 

W = ~w 1 r 4 fg, (2-8) 
(vro, veo) = (l-^r, , (2-9) 
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l 2 = \r 2 g, (2-10) 

where n\ and W\ are the number density and enthalpy in front of the shock. Since 
we consider a strong shock, the external medium is cold and non-relativistic, and 
thus the pressure of the medium in front of the shock can be neglected. We also 
assume that the density n\ varies as a power of the radius: 

ni cxr- k . (2-11) 

The functions /(%), g(x) an d h(x) can be obtained by solving differential equations, 
but they reduce to simple analytic forms for the special case m = 3 — k as 

/ = x -(4fc-ir)/(3fc-i2) j (2 . 12a) 
9 = X~\ (2-12b) 

fc = x -(2*-7)/(fc-4) i (2 . 12c) 

We only examine the stability of the case that includes the point explosion into the 
constant density k = 0. The calculation is simply, and the result is likely to be 
common to other cases, in which the growth/decay rate may be modified. 
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2.2. Perturbation equations 

We denote the Eulerian variation of a variable by S, and we expand variables, 
e.g. the pressure as 

p = Po + 5p = Po {l + (r 2 ) s Z P P l } , (2-13) 

where Pi(cos9) is the Legendre polynomial of order I, and a power law form is 
assumed with respect to the "time" r. Since the background quantities change in 
a self-similar way, the perturbations also have self-similar forms. The index s is 
a complex number in general. The real part of it represents the growth rate for 
5ft(s) < 0, since the factor r decreases with time. In a similar way, perturbations of 
the density and velocity are expanded as 

D = D + 5D =D {l + (r 2 ) s ( D P l } , (2-14) 

Vr = v r0 + sv r = i-^{i + (r 2 y^ R p l } , (2-i5) 

v e = Sv eo = (r 2 r^T-^Pi • (2-16) 

A perturbation of the function W can be expressed by the perturbations of the 
pressure and radial velocity as 

W = W + 6W = W {l + (r 2 ) s g(2^ P + ^ R )P l } . (2-17) 

Using the definitions in Eqs. (2- la) - (2- Id), we have the linearized perturbation 
equations. After separating the time and angular parts, we have the differential 
equations for the functions £d,£p,£r and £t- The explicit forms can be written as 



(k - 4) X ^ - \{k - 4) X ^ = (k — 3Ud + \{k- 2)U 

-i{-3(s + 4) + fc( S + 6)Kp, (2-18a) 
(k-A) X ^ + 3 -(k-4) X ^ = -\3s + 5 --k(s + l)}^ 



d X 2" //x d X I 3 V 3 

+ I {3( s - 4) - k(s - 6)}£p, (248b) 

(*-4)xf = -{ 3a + ^- fc ( a + I)} fr> (2.18c) 

(k - 4) X ^ + 2(k - 4) X ^ = (k - 3U D - 2(3k - 7)& 

-1(1 + 1)Zt, (2-18d) 

where we have retained the highest order of r only. The differential Equations (2- 18a) 
- (2-18d) can be solved analytic ally. The solutions are expressed with four integra- 
tion constants, a, b, c and d: 
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= a X Pl + b Ql I (I + 1) X P2 + cq 4 _ X p+ + ^4 + X p - 
Cfl = bq 2 I {I + 1)X P2 + «/ 5+ X P+ + ^5_X P - , 
Z P = bq 3 l(l + l)x P2 +cx P+ +d X p - , 

= b X p ' 2 ■ 

Here the indices pi and coefficients % depend on s only: 

fc-3 



Pi 

P2 
P± 



k-A S ' 

A;(3s + 7) - 9s - 14 

3(jfe - 4) ' 
-fc(6s + 17) + 18s + 43 ± /i 
6(fc-4) 



(2-19a) 
(2-19b) 
(2-19c) 
(2-19d) 



(2-20a) 
(2-20b) 
(2-20c) 



where 



and 



5i = - 



/i = {fc 2 (48 s 2 + 192 s + 84l) - 2 fc (l44 s 2 + 528 s + 1943) 
+432s 2 + 1440s + 4489} 1/2 

7 (As — 2) [2k 2 (4s 2 + 23s + 5) - k (48 s 2 + 240 s + 43) 



(2-21) 



+72 s z 



(2-22a) 



<?2 



306 s + 46}] {k 2 (4 s - 2) - k (26 s - 5) + 42 s - 2} 

= -{2fc 2 (4s 2 + 23s + 5) -fc (48s 2 + 240s + 43) + 72 s 2 + 306 s + 46}"* 

x {fc(4s + 1) - 2(6s + 1)}, (2-22b) 

^3 = {2 k 2 (4 s 2 + 23 s + 5) - k (48 s 2 + 240 s + 43) + 72 s 2 + 306 s + 46}"* 

x{2(2/c-3)}, (2-22c) 

^ 4± = [4{fe (s + 4)-3s-10}{fe 2 (4s 2 + 9s-23) - k (24 s 2 + 50 s - lOl) 

+36 s 2 + 69 s - lio}] ~* [k 3 (l2 s 3 + 113 s 2 + 231 s - 406) 

-k 2 (l08 s 3 + 949 s 2 + 1822 s - 2736) + k (324 s 3 + 2643s 2 + 4763s - 6126) 

-324 s 3 - 2439 s 2 - 4128 s + 4556 ± [k 2 (s 2 + 3 s - 14) 

—k (6 s 2 + 17 s - 62) + 9 s 2 + 24 s - 68} /1] , (2-22d) 

g5± = [8 { fe (s + 4) - 3 s - 10}] _1 (-29A; + 67 ± /1) . (2-22e) 

The physical interpretation of the solutions is evident. The term with the constant 
a simply represents the perturbation of the number density, without disturbing any 
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other quantities. The terms with the constant b represent the perturbation of the 
angular momentum. This mode is possible only for non-spherical perturbations 
(Z ^ 0). For a spherically symmetric perturbation, £r = 0, Z = 0. In this case, the 
solution should be reduced by setting b = in Eqs. (2- 18a) - (2T8d). The terms 
involving the constants c and d represent the pressure waves propagating in the 
inward and outward directions. 

2.3. Regularity condition at the origin 

For the boundary condition at the origin, we require that the fluid does not 
undergo divergent perturbations. Also, the pressure perturbation should vanish: 
bp — > . This condition corresponds to 

£p -> (2-23) 

for x — * 00 • It restricts the possible form of the power index K(s); that is, £p must 
decrease as the origin is approached. The explicit condition is written as 

2.4. Boundary conditions at the shock front 

Here we consider the deformation of an expanding shock wave. The radius of 
the shock is assumed to be given by 

r = *i 1 -57 — rr^l + f. ( 2 ' 25 ) 



2(m + l)r 2 

where the first term denotes the radius of the unperturbed shock and the second 
denotes the perturbation. We also assume that the perturbation function has the 
power-law form 

CY / o\ s ~ 1 

v = --t[r 2 ) P l} (2-26) 

where a is the normalization factor. This change of the shock front induces a radial 
velocity as 

W =§ = - a i^(r^-'p,. (2.27) 

We now apply the jump condition for the perturbation functions at this perturbed 
shock front. That is, we stipulate that the number flux, energy flux and momentum 
flux be continuous across the surface (see |A|). The values at the perturbed posi- 
tion are evaluated by the Lagrange perturbation. For example, the pressure at the 
perturbed shock position is given by 

Ap = 5p + r?^ . (2-28) 
or 

It is also easy to rewrite the boundary conditions at shock front as 

£ D /a = 21s - 10k + 13, (2-29a) 
£ P /a = 9s-4k + 5, (2-29b) 
in/ a = 12s -bk + b. (2-29c) 
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2.5. Approximate solution 
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Fig. 1. The growth rate as a function of normalized tangential velocity, |{t/«| • l{l + 1). Two curves 
are shown corresponding to the external density distribution k. Uniform case is k = 0. Dots 
near K(s) = —1.5 are the maximum values which are determined from the regularity condition 
at the origin. 

We can easily determine the integration constants by matching the jump condi- 
tions to the interior solution Eqs. (2- 19a) - (2-19d) (see||). The temporal dependence 
s is specified by one parameter, the velocity in the tangential direction. We use the 
normalized value |£t/o| -1(1 + 1) and plot the growth rate 3ff(s) as a function of it 
in Fig. 1. Both the cases £t/o < and £t/« > are shown. For the spherically 
symmetric case £t/& ■ 1(1 + 1) = 0, the growth rate is 3?(s) ~ —12.7 for k = 0. As 
the value £t/« decreases, 5ft(s) increases. We only show the region of 3?(s) limited 
by Eq. ( gjg ); that is, the curves are terminated near 3?(s) = —1.5. As the value 
£r/a increases from zero, 3?(s) monotonically decreases. This implies that the tem- 
poral behavior depends significantly on the r factor. The mode grows strongly with 
the deceleration of the background fluid. In particular, the growth rate is large for 
shorter wavelengths in the tangential direction. In this figure, we also show the effect 
of the external density gradient, uniform (A; = 0) and a decreasing slope with the 
radius (k = 1). The general features of the growth rate do not depend on the density 
distribution. The growth rate is modified by about 10%. 

The radial functions of the background and perturbation fields are shown in 
Fig. 2. We display the density profiles near the shock front, which corresponds to 
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1 

,. 2. The shape of background density(Do) and its perturbation value((5_D) normalized at shock 
front. The parameters adopted here are Bl(l + 1) = 1000 and k = 0. \ is a similarity variable 
which show the distance from shock front, where x = 1 is the shock front. 
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X = 1. The thickness of the expanding shell is approximated by the region from x = 
1 to 2. The radial function of the perturbation is localized within the small region 
of the shell; that is, it is a quite steep function. 

§3. Summary and discussion 

In this paper, we have studied the linear stability of a spherically expanding 
shock in the ultra-relativistic limit, i.e. the large r limit. The temporal behavior 
of the perturbation is expressed by a power law in r. We found a solution that 
grows with the deceleration of the shock. The growth rate becomes large for short 
wavelengths of the angular fluctuation. That is, there exists a rapidly growing mode 
for the non-spherically symmetric displacement. VishniaclllP discussed the physical 
interpretation of the unstable mode in the non-relativistic blast wave. It is important 
to remember that two kinds of pressures exerted to the shock plane have completely 
different natures. In the exterior region before the shock, the ram pressure domi- 
nates and is parallel to the direction of the shock propagation. After the shock, the 
thermal pressure dominates and is isotropic. When the shock spherically expands, 
two directions, i.e. the direction of the propagation and the normal direction of the 
plane coincide. However, when the shock is rippled with the tangential velocity, two 
pressures lose the balance in the tangential direction. This unbalance is expected to 
be larger with the increase of the velocity. The unstable mechanism indeed works 
in the ultra relativistic blast wave as shown in this paper. This result depends on 
several assumptions and approximations, but it is quite suggestive. Suppose, e.g., 
for the relativistic flow dynamics from r = 10 2 to 1 associated with the gamma-ray 
bursts, that the amplitude increases by more than a factor of 10 20 . The dominant 
scale of the unstable mode depends on both the initial disturbance and the growth 
rate. If there is no characteristic wavelength in the initial stage of the expansion, the 
unstable mode is likely to be non-linear at the smaller scale due to the large growth 
rate. This mode may lead to a number of fragmentations. The radiation from the 
blobs may be reflected by an irregular time-profile of the gamma-rays. 

In order to check this possibility, further study is necessary. For dynamical 
evolution, we must simulate relativistic flows using axi-symmetric or 3-dimensional 
numerical codes. Fluctuations are manifest only in the small region after the shock 
front. The structure may be found only by using a fine resolution in the numerical 
calculations. The unstable mode may be suppressed to a certain value of the am- 
plitude. Also, the size and number relations of the fragmentation can be estimated. 
The effect of the flow dynamics will appear in the emission, but the gamma-ray emis- 
sion process is not clear at moment. However, the complex structure of the shock 
front may universally give the time history of the burst, irrespective of the detailed 
behavior of the emission process. 
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Appendix A 

Jump Conditions 

The difference between the values on the two sides of a surface is denoted by 
[a] = a\ — 0,2, where the subscript "1" indicates the fluid located before the shock 
wave, and "2" after the shock wave. The number density and energy-momentum 
are continuous across the shock surface. The jump conditions for the relativistic 
adiabatic shock are given by 

[n/JV s J=0, (A-l) 
[T^N s/ ,N sv ]=0 , (A-2) 
[T» v N 8ll U av ]=Q, (A-3) 

where T tIU and are the energy- momentum tensor and the 4- velocity of the fluid. 
The time-like unit vector N s ^ describes the motion of the shock front, and the space- 
like unit vector U s ^ describes the normal direction of the surface. For a spherically 
expanding shock wave, they are given by N SfM = {rV,T,0} and U SIX = {r, rV,0}. 
We consider linear perturbations of them. The explicit forms are given by 

N Sfl = {r (v + r 2 5v r ),r (1 + vr 2 5v r ) ,n ± } , (A-4) 
u 3(1 = {r (i + vr 2 5v r ) ,r (y + r 2 5v r ) ,u ± } , (A-5) 

where 5V r is the velocity perturbation in the radial direction, and A^j_ and U± are 
the tangential components of the perturbation of the variables. Their explicit forms 
are not necessary in the calculations. 

Appendix B 

Analytic Solution 



The coefficients of the analytic solutions Eqs. (2- 19a) - (2-19d) are determined 



as 

a = [28 (k - 2) jfc (4 s - 5) - 12 s + 7} 



x 
x 



\k 2 (4 s 2 + 9 s - 23) - k (24 s 2 + 50 s - lOl) + 36 s 2 + 69 s - lio}] 1 
-8 k 5 (352 s 3 - 584 s 2 - 3983 s + 5905) 
+k A (5616 s 4 + 22872 s 3 - 115057 s 2 - 225180 s + 425829) 
-k 3 (60624 s 4 + 56344 s 3 - 794477s 2 - 518881 s + 1515622) 
+k 2 (242352 s 4 + 11088 s 3 - 2362877 s 2 - 277212 s + 2660875) 
-k (423792 s 4 - 125064 s 3 - 3211707s 2 + 457097s + 2302554) 
+2 (136080 s 4 - 61020 s 3 - 817893 s 2 + 227748 s + 392500) 
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-3 {8 k A (6 s 2 -s- 15) - k 3 (l08 s 3 + 451 s 2 - 420 s - 867 
+ k 2 (900 s 3 + 1418 s 2 - 2409 s - 2273) 
- k (2484 s 3 + 1605 s 2 - 4737 s - 2558) 
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+2268 s 3 + 342 s 2 - 3060 s 



1040} h 



(B-la) 



4 Z(Z + 1){A; (4s - 5) - 12s + 7} 



-1 



8/c 3 (20s 2 + 155s-38 



k 2 (384 s 3 + 4012 s 2 + 7249 s - 1265 



+ k (2304 s 3 + 17844 s 2 + 13577 s - 1539) 

- 18 (l92 s 3 + 1208 s 2 + 439 s - 30) - {8 k 2 (2 s + 3) 
—k (36 s 2 + 125 s + 63) +2(54s 2 + 75s + 20)}/i] , 

2 {7 - 12s + k (-5 + 4s)} {2k 2 (4s 2 + 23s + 5) 
—k (48 s 2 + 240 s + 43) + 72 s 2 + 306 s + 46}] ~ 
x [{8 k 2 (2 s + 3) - k (36 s 2 + 125 s + 63) + 108 s 2 + 150 s + 40} 
x{-2/s 2 (8s 2 + 34s-2l) + k (96 s 2 + 372 s - 149^ 

- 144 s 2 - 504 s + 121+(2 k - 3)/i}l , 



(B-lb) 



d = 0. 



(B-lc) 
(B-ld) 



[i] 

[2] 

[3] 

[4] 

[5] 
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